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ABSTRACT 


We have used Monte Carlo simulation, lattice dynamics in the harmonic 

approximation, and solution of the hypernetted chain equation to study the 

classical two-dimensional one component plasma. We find a fluid phase for 
— . 2 i 

L = e (rm) AgT < 125 + 15 and a solid phase for higher T. The solid phase 
shows directional long range order. In the solid phase positional long range 
order is lost as the thermodynamic limit is approached. We also present the 
results of calculation of the thermodynamic functions and one and two particle 


correlation functions. 
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This paper is concerned with the properties of a two-diinensional one 

I 

component plasma. Our system consists of a single species of charged par- 
ticles immersed in a uniform neutralizing background. The particles interact 
via a 1/r potential, where r is the two-dimensional separation. Our calcula- 
tions are limited to ranges of temperature and density such that quantum 
effects are unimportant. We have made calculations of the equation of state 
in the' fluid phase using both the hypernetted chain equation (HNC) and 
Monte Carlo simulations. Our calculations in the crystal phase were done 
by Monte Carlo methods. 

There are two reasons why we find such a system interesting. First 
it can be considered as an idealized model of a bound surface layer of elec- 
trons above liquid helium four. Second there have been extensive simulations 
of the properties of the three-dimensional one component plasma. The exten- 
sion to two dimensions may provide insight into the behavior of both systems. 

We begin by briefly reviewing the state of our knowledge of the electron 
surface layer above liquid helium four. Several years ago Crandall and 
Williams^ suggested that under favorable circumstances electrons trapped on 
the surface of liquid helium might crystallize to form a two-dimensional 
electron solid. Since in most experimental situations the density of electrons 

can be changed by several orders of magnitude, it was hoped that the so called 
2 

Wigner crystal might be within experimental reach. This led to a great deal 

3 

of theoretical and experimental activity in the following years, and Chaplik 

suggested that a similar crystallization can occur in the inversion layer 

near the surface of a semiconductor. In the helium context a model of a charge 

compensated one component system of N electrons confined to an area A at a 

temperature T interacting with 1/r potential can be and has been considered 

4-7 

as the canonical model. 



, In this paper we consider only the classical behavior of the model. 

This is appropriate for the electron surface layer above liquid helium 
in the usual experimental regime. However, our model is not appropriate for 
the problem of the metal -oxide-silicon inversion layer where the electrons 
form a degenerate quantum system. 

3 

.41though studies by Brown and Grimes of cyclotron resonance 

in a tipped magnetic field have shown that the electron motion 

on the surface of liquid helium is two dimensional, it is clear that 

for a strong clamping field (most experiments require a clamping field in order 

to localize the electrons layer for a reasonable amount of time) one needs in principl 

to take into account the coupling in the perpendicular direction, for example 

9 

the deformation of the helium surface. However, the characteristic dimensions 
are such that the interelectronic spacing ( ~ 10 A) is much larger than the 

9 o 

spread in the charge density in the direction perpendicular to the surface (^10 ■“A) 
so that the system can be considered to be essentially two-dimensional. There- 
fore the model of a two-dimensional electron gas, neutralized by a uniform 

2 , 

positive background, and interacting by a e /r potential is probably a reasonable 

first approximation to the experimental situation. The system is characterized 

2 4 

by the dimensionless quantity e /ak T, where a = (irn)^ and n = N/A. 

The simulations that have been made on the three-dimensional one component 
plasma have established its equation of state, the phase boundary between the 
crystal and liquid phases, and the two particle correlation function 
at several densities and temperatures. In addition, Lindemann's ratio at 
melting has been found to be 0.17— — rather close to the values for other inverse 
power potentials. The height of the first peak of the structure factor at 

freezing was found to be close to 2.85 again close to the values for other 

potentials. In view of these results it seemed worthwhile to carry out a 
similar study of the two-dimensional one component plasma. In particular we 
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were interested in determining whether the two-dimensional system would 

undergo a phase transformation to a crystalline phase. To our knowledge 

there has not been any calculation comparing the free energies of the solid 

and the liquid phases which is* after all, fhe basic method to locate this 

phase transition. In this paper we present such a calculation. We employ 

both the hypemetted chain integral equation and the Monte Carlo technique 

to calculate the free energies. We have computed the thermodynamic functions 

and correlation functions over a wide range of T, 1 < T < 300. A recent 

publication by Totsuji^^ contains Monte Carlo results for the thermodynamic 

functions and pair correlation function for 0.15 < F < 50. Within the 

range of F our results are generally in good agreement with those of Totsuji. 

On the other hand our results are quite different from a very recent computer 

experiment^^ which employed a special type of molecular dynamics method 

(PPPM: Particle-Particle/particle-Mesh) . Contrary to the X point transition 

obtained there we tentatively find a first order transition, our transition 

being roughly 20% higher in F; namely F = 125. Our results are qualitatively 

12 

similar to the corresponding three-dimen.sional calculations of Hansen and 

13 

Pollock and Hansen . We find that the triangular lattice is stable and have 

calculated the harmonic phonon dispersion laWs for such a lattice. We also 

find that the two-dimensional square lattice is dynamically unstable. Our 

calculations for the harmonic solid are in agreement with a recent calculation 

14 

by Bonsall and Maradudin 

Before proceeding further we should comment on the -existence of two- 

15 

dimensional crystalline order. Some years ago Mermin published a rigorous 

proof, based on Bogolyubov's inequality, that two-dimensional systems cannot 

display long range crystalline order. The proof had two limitations. First, 

and probably less important, the interaction po-tential was assumed to fall 

2 

off faster than 1/r . Second, the result only applies in the thermodynamic 
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limit. liVhen the same mathematical methods are applied to a large but finite 
system one finds that no inconsistency arises from the assumption of crystalline 
order. Thus any system that can be studied in the laboratory or in a computer 
simulation can exhibit crystalline order. We do indeed find a stable crystal 
phase in our simulations, however, it does have unusual properties. These 
are described in sections VI and VII. 

The plan of our paper is as follows, after formulating the problem in 
section II, sections III and IV are devoted to the calculations of the thermo- 
dynamic functions in the liquid and solid phases respectively. In the liquid 
phase we present results for both the HNC method and MC simulations. Section V 
is devoted to the determination of the phase boundary between the crystal 
and fluid phases. Lindemann's ratio and its dependence on the size of the 
system are discussed in section VI. The one and two particle distribution 
functions are presented in section VII. 


II. FOBMULATION OF THE PROBLIM 

Consider a system of N electrons obeying classical statistics, confined 
to an area A, and neutralized by a uniform positive background. The Hamilto- 
nian, apart from the kinetic energy, which does not enter into our considerations, 
is then 

2 

H = S -2- - jNnv(O) (2.1) 

i< j ^iJ 

where, the last terra, arises from the interaction of electrons with the uniform 

I ^ I 

positive background. Here r, , = r. - r. is the two-dimensional separation 

ij ' i j ' 

9 

and n = N/A is the area density of electrons. The Fourier transform of e~‘/r 
potential in two dimensions is given by 
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v(k) a 8 j2-_ 4^r » -r— 


( 2 . 2 ) 


Once again k a (k ,k ) la a two-dimensional vector. For the problem of elec- 
X y 

trona bound to the autface of helium four one should use the *' /enomalitjed'* 

7 

Charge 


e ^ “ efface, + 


<2.3) 


where we have assumed that uniform media with dielectric constants and e., 

X. 

fill the adjacent half spaces. The thei'modynamics of this one component 


classical electron plasma is determined by a single dimensionless parameter, 

^ 2 

r « rV- . 

k^Ta ’ 

B 


Where k la the Boltzmann conatamt, T the temperature and a ia the average 
B 

2 

interparticle distance determined by a « l/(na ) , The complete solution of 
the problem therefore reduces to the calculation of averages of the kind 


. jV(r^...rj^ekp t - S v(r )) d^r^. . .d^r^^'Q 


B‘i<j 




(2,3) 


xt?here 


“ J • * ‘J i If T 

"* i<j 


fir • .2 

^ ver, ,) f d r. . . .d . 


IJ' 




The two particle distribution function is given by 


P(r^,r 



Q« ^ 


.jexp 



- v(r d^r^. 
Kj 



(2-.S) 


^ ^ , 2 

For the liquid phase » n g(r^ 2 ^> where g(r) is the pair coi'X-elation 

function. 
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Before we proceed any further we should point out that the stunmation of 

the ring diagrams (random phase approximation) produces a divergent free energy 

and therefore the result analogous to Debye-Hiickel limit in three dimensions 

does not exist. In fact it is easy to mimic the corresponding three-dimensional 
IS 17 

calculation ’ to obtain 


^ = N log(l/n) + 1 + W 
r*i • 
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(2.7) 


where W is the contribution of the ring diagrams 


W = W 


2D 

ring 


1 

4TT n 



r 2rme 
\ k^T 


2 


i - iog(: 


/'2TTTie^ 



( 2 . 8 ) 


which diverges as -* «. This is different from the behavior of the corres- 
ponding three-dimensional integral; 



This high q (small r) divergence for the contribution of the ring diagrams 

18 

to the energy has also been noticed by Totsuji . He has, however, shown 
that if one includes the simplest set of next order diagrams then the diver- 
gence is cancelled. His results can be expressed as a small F expansion 
for the excess internal energy, 


ex 


NkaT 


= T + 2y - 1 + 2070) + 


r^C-8(27T)^ + 8(1 - 2 y)M" + 4(1 - 2y) 2r2] -f 


(2.9) 


where y is Euler’s constant. 



Ill, LXQtJID PIUjSE CALCUl^VX'XOjrS 


aiveii the paii' coi'relafelou function fix*) ehcn the excels internal, enerfy t| 
ia fiven by, 


>» 




.ex 


Nk r ' “ a J dkisik) - 1] 

‘Bo o 


l3»l> 


where Blk) is the Structure factor defined by 


«» 


s(k) « 1 a 


rd V kr ) Cf i x*) ‘*'11 dr 
o 


\3,a) 


where d xx) Is the zeroth order Bessel functioui, the pressure obtained from 
o 

tire ririal theorem is 


pA 


* 1 4 * 


U 


ex 


k^T 2Nk^t 


13.1) 


and excess free energy is giren by 


y pro 

^ dC 

*36 




o 


13.4) 


the ideal gas fx*ee energy is given by 




, .V - 1 


j?k^T \m^ ^ 


13.5) 


or , in ovxr units , 


k t X 

B 'V 


13 . 6 ) 


4 


where * e in, ah“ is the three-diaienaional Rydberg’s oonatant, 


3 


Equations (3.1) to (3.8) are the relations we use in our liquid phase calcula- 
tions . 

A* The Itypernetted Chain Method 

The hypeitvected chain equation has been very successful in predicting 
the properties of the danse one component plasma in three dimensions, in 
this subsection we shall present the results for the two-dimensional case. 

The hypernetted chain integral equation for the pair correlation function is 
defined by the following two equations 


g(r) - 1 c(r) + njdr' Cg(r') - l]c(jr - r' 1) 


( v) W ) 


^ g(y) - g(t) - 1 ~ c(r) - 


v(r) 

V 


(3.3) 


where c(r) is the direct correlation function. This equation is difficult to 
solve for long range interactions, Following the method used recently by 
Springer, et al“ we carry out a subtraction procedure to rewrite the equations 
in terms of short range quantities. W& decompose the original potential 
v(r) into long and short range parts 


v(r) ^ V- (r) + V (r) « ei*f(br) + - (1 - Qi*f(br)] . (3,9) 

ir sr r r 

The parameter b is for the moment arbitrary. We denote the sum of nodal 
diagrams by *N;(x*) ; 


- g(t) - 1 - c(r) 


(3.10) 
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and decompose N(r) into short and long range parts 


Ngj.(r) = N(r) - v^^(r) 


( 3 . 11 ) 


Since at large r, N(r) F/r, is short ranged. Similarly for the 

direct correlation function c(r) we have 


c (r) = c(r) + V (r) . 

sr ir 


(3.12) 


It is now easy to see that Eqs . (3.5) and (3.6) can be solved by Fourier 
transforms. The set of equations to be solved are 


c (k) [c (k)-^Erfc(k/2b)3 - ~Erfc(k/2b) 


sr 


sr sr 


1 - [c (k) - -^Erfc(k/2b)] 
sr r 


(3.13) 


g(r) =s expCNg^(r) - 


(3.14) 


c (r) = g(r) - 1 - N (r) 
sr sr 


( 3 . 15 ) 


where 


^ W * T* ^ 

c (k) = n I e c (r)d“r = 4ir ! J (kr)rc (r) dr , 
ar J sr J o sr 


(3.16) 


and 


N (k) 
sr 


p ik*r N (r)d“'r = Irri J (kr)rN (r) dr 
= n e sr <i o sr 


(3.17) 


We also have for the inverse transform 
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s s r '* V 

o 


kN (k>J (kr) dk . 
o 


Xt: is vjlear Sroni Sds* (3,T - 3,16) thar the original equation has been 
cast intc a i^orav which involves onlj' short rangre functions, and the numerical 
traasSorms can therei'ore he calculated much more readil>% The method is 
siailar to the E\vald technique for calculatinif lattice sums, and depends on 
an optimuw choice of the parametex' h such that the functions are short ranged 
in hoth k and r space* The remaining problem is to calculate the Hankel 
traxis,foj'ms accurately. Since we will need to take transforms hack and forth 
hetween r and k space, it is important that the algoritJim preserves the ortho- 
gonality of Bessel functions. Such an algoritlua has been worked out hy 
'"'X 

Lado . A significant impj’ovement of the accuj'acy was obtained by calculating 
the iiaros of the Bessel function to machine accuracy. Machine memory lia^ita-^ 
tions esaexxtially restricted us to a maximum of 200 points in r and k space, 
Oux* x’esults are presented at the end of this section, 

B. Monte Cax'lo Method 

Oux* MC simulations were perfoitned in the time-honored mannex* pioneered 
hy Metropolis et al'*'", Eacix run had a given number of particles N, a given 
area A* and a fixed value of T. The area A was chosen to be a rectangle 
capable of accomodating a section of triangular close-packed lattice with 
periodic boundax*y conditions, To minimise surface effects the rectangle was 
chosen to he neax*ly square with the ratio of the x and y edges as 
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Due to the long: range nature of the potential, the interaction of each 
particle with the other (.N - 1) particles in the basic rectangle, with all 
images of the particles, and with tlie uniform neutralising background must 
all be included. This complete interaction can be written as 


'•(?) = f 7= — =-Trr 


L h 

X j 


d" P 
r -t p 


(3.20) 


where X* as the set of vectors generated by the basis vectors (1,0) and 


(0, y'of'2) and where r is the distance vector between the two particles. This 

23 

interparticle potential may be efficiently handled via the Ewald techniq.ue~ . 

1*> 

The computation proceeds in a similar way to the 3D case “ . The result is 


erfc(va!r.L 4. X‘ |) - „ 

v(r) » X ~ 

Vx' 


X exp(^2TTi 


r/L + X^ 
X"-? , 1 




X” \r 


1 . -ir’L 

ertc(-- 


(3,21) 


where the second sum is over a reciprocal lattice with basis (1,0) and (0,2,v 3) 
The prime on that sum denotes the exclusion of the term with .V' = (0,0). 

The parameter a may be varied to achieve a balance in the rates of convergence 
of the two sums . 


Even with the application of the Swald technique equation (3.21) is 
unacceptably slow for .MC calculations. To make the calculation more efficient 
v(r) is split into two parts. The first part consists of the spherically 
s^Tiunetfic (X* =» 0) term which is tabulated with a 35,000 point mesh. The 
remaining part of v(r), which is invariant under reflections (x “X, v -y) 

and inversions (r -* -r) , is tabulated on a 151 x 1”1 point grid. With linear 
interpolation applied to eacli part of v(r) the potential energy of a configura*- 




if*, 
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tion may be efficiently calculated with an error of approximately 0.001% , 
which we found to be negligible. 

If we now include the interacxion of each particle with its own images 
we obtain the total interaction energy of a configuration of the X particles 
and images as 



v(r.j) 



■ rivuiiNA-ly i 
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(3. 


2 2 ) 


The constant term is just half the Madelung energy of a rectangular lattice 

with sides L and L . In our units 
X y 



-1,09653 r f — 

\L L / 
X y 


(3.23) 


In our MC calculation the center of mass is not fixed. To correct for 

this effect the difference between the MC value for the excess internal energy 
MC 

IT /NkT and the static energy U^/NkT of a perfect triangular close-packed 

24 12 

lattice must be multiplied by N/(N - 1) . If we follow Hansen and call 

the difference between the excess internal energy and the static energy the 
thermal contribution (U ) to the excess internal energy 


_ fD-'"' “on N 

NkjT M-i 


Our excess internal energy is then 


(3.24) 


ex 




U 

NICbT 


th U 


NICbT 


(3.25) 



Here the atatic lattice energy ia given by Unli the Madelung energy oi the 
triangular oloee~packe€l lattice. 

n -1.106103 r . (3.38) 

Nk„T 

B 

The liquid phase Monte Ciurlo runs were made with 16, 36, 64, and 100 
particles. These numbers of particles, while small for three-dimensional 
work, are reasonable for two dimensions. The starting configurations for 
the runs were either a triangular close-packed lattice or a quasi-random 
configuration. For each run approjvimately 13,000 moves per particle were 
made and on the average 50% of these moves were accepted. 01 these approxi- 
mately 3,000 movea per particle were disoai'ded in order to allow the system 
to lose its memory of its origihal configuration and reach an "equilibrium 
state." Since the amount of cotifiguration space to be sampled is considerably 
reduced in two dimensions over three dimensions these runs represent very 
long Markov chains when judged by standards normally applied la three dimensions 
Huns of this length are needed to accurately determine the thei*mal portion of 
the excess Internal anex-gy which is Only about 1% of the total excess internal 
enex'gy. From block aVQX'aging we estimate that our calculated excess thermal 
energies have errors of the order of 1%. 

For the amall F region effects due to the small numbers of particles are 
not really significant; even the 16 particle rxina show little effect due to 
the small siKe. For the largest values of F number dependence becomes much 
more impox’tant. However, as can be seen from Fig. 1 the number dependenca, 
in the energy, is essentially eliminated by the time 100 particles is reached. 

The most convenient way of detxling with the HKC and MC results is 
with a simple and accunxta fitting function. 



u 
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Ij 
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Tha valvift ajg u, and a.U othar ilttinir payamttara, ia fivan in Sabla I. Wa 

bava a aaoond iittin^* tojnnuia iot valnaa oi ^ in t;ba i'nn«a 0,45 to 130. It 

12 

la baaad on Ttanaan'a tb5*aa~d4maniional woi‘K and tka fact that 


Li»\ 






St «1.1Q6 f* 


(3,27) 


Tha f Omni a la 




(a I !■ 


T’>\ 



n’ 



^7 "I 


(3.2$) 


whai^a tht value o,t ia lifted by aquation (3,2$), Thia fiotsnula acourataly 
pepx'OdUQ^a both oux’ MC and MC tesulta with the pax'ametarp jsiven in table I. 
Ouy MC teanltp ana ps-eeanted in table 11 . and 5i§:ui'a 1, 

We have oalouiated the a^toeae fx'ee enevfy o.t the liquid phaae hy uaing 
equation (3,4) with equations (2,9) and (3,2$), We awltoh between equation 
(2,9) and equation (3,2$) at P^ « 0,43, Ou,n li'ee enets)’* is then siven by 


ax ,.,p « 

m ^ X -h M 3 

Vs '». Lh ^ 


T 


^ (3«ayw + (i-av)(j^a-i) « | 


(3,29) 


oP^ 


tox‘ 0 < r < 


and hy 
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+ a 


+ a 


f!i 

f . 1 - - 


- f - ^ ^ 

iL 2 

\ 2 

\ag+n- 


\l+a F F /J 

o 6 o 


f - :A - 

1 ' 


’L 3 





(3.30) 


tor r < r < 130, where T is given by equation (3.2 9) evaluated 

O'”—* O -B 

ior r = r^* Our results for the free energy are tabulated in table 2 and 

displayed graphically in figure 2. We agree well with the results of 
10 

Totsuji for the free energy. It should be noted 

that our Monte Carlo runs are approximately ten times longer than his. For 
this reason we believe our data ia more accurate. 

The excess specific heat at constant area may be calculated by differen- 
tiating the excess internal energy 


% , fi 


B 


ar LF NkgT. 


(3.31) 


In the region 0 < F < 0.45 we obtain 


,ex 


Nk 


'B 


= - F^2Y 


■f 1 + Sti 4 + 2 S/n 


■m 


+ F^C(2Y”1)4(2+3 ^ 2) + a(6Y-l) ^ F + 24(^ F)^]-5cF® 


(3.3 2) 


and for 0.45 < F < 130 we get 
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S', r^_!i_ + ^”3 

’”'b -a^ + r (aj + (a,, + (a^ ♦ D'‘ 


R n 

5 , r 1 

■*t* 


ra2r ^1 ^ 

k ■ — — -— + 


(3,33) 


_ 1 


r (<i_^ D (a^ + n'" (a i- tY 


The equation of atate follows from the virial theorem, equation (3.3). 

Table 2 contains out results foi' the specific heat* Fi^ux*e 3 shows our 
results for G_^ usin^ the flttinf formula and from direct Monte Carlo calcula- 
tions. 'The latter calculation is difficult due to the large fluctuations in 
the data and our I’esults are preliminary in nature. We will discuss them in 
more detail in section V. 

IV. SObJD PHASE CAICULATIQNS 

A combination of lattice dynamics and MC methods was used to compute the 
properties of the solid phase. The lattice dynan^ics calculation was fjerformed 
in the manner of Bonsall and Maradudin , The technique of special points" 
was used to efficiently calculate various chermodjuiamic functions by averaging 
functions of the frequencies over the first Brillouin sono. The lattice dynamics 
calculations provided us with approximations for the free energy and ocher 
thermodynamic functions and the root mean square deviation, of particles froia 


^quHibx’luttv Xuttlcse Tb«? MC ci^ioulationa pvovivied ua with 

ibtex'mX a»d the om* anti two pantioXe 4iatyibution funotiona 

fox' the ox'jratalXine phaae* ftow these x'eavxlta we oalouiateu the fx'oe energy 
and ape«l.tlo heat. These MC oaXonlatlona we.te pei'tox'mevi in the aame way 
aa in the liduivi phase with the exoeptlon that All rvma wei-e atax'teO fx't'iw a 
pex'teot iattioe, Again nma wex'© pex'tovnved for lA 3d, dA, aixci 100 pax'tiolea. 
fox' r < ISO we tovxxxd that the deviationa ot the pax'tivJiea fx'Ottv theix' origlxxal 
lattloe aitea divi not i«each a steady yalne vhu'tng n,.na. This value o£ was 
taisesx as a px'eliwinary Indioatlon of the looatioxx of the phase ti'aneition to 
the liv\vxid phase. Solivi phase MC x'uns were saatie for values of r between 130 
axni dOO. 

A oeixtex' of mass ox'vreotlvHV was applievl just as Ux the Hviuid phase, the 
depexxdenoe of the tUermovivnamlo v\uantities vvn the nuiwber of pax'tioles In the 
s.tftxiilatioxx was foiuvd to be less sevex’e thaxx in the high T liquid phase with the 
100 pax'ticle syatew again well X'epx'esenting the inf ixxite pa,x*tiole system. The 
pyediotlon of harnvonic lattice dynaittlcs for the excess Ixxternal energy is 


^|Harx« 




id.l) 


the v>X'essure is, of coux'Se, still given by eu\xation vd,.' i. 

lx\ the harmoniv' apprtui4t\ativ>n the total free energy is given by 


pHarw 



Tnwtieiiiii je 


r ‘Ib^Vd'' 

sixth 

«N>i 




whex'o the Sixm is vwer fx'exiuexnxies tn the fix'St x-oixe and j is suwxnxed over the 
two bx'axxchea. In the olsssioal liaxtt we get 


^harro 



-A 



v-n X 


H,, 

w vV 
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In our units we obtain; 


OBIGWAL » 

OP POOB 




arm U 


Nk T Nk_T N 
B B qj 


OD (q) 

liJ. 


_ -1 w . V 'ly ^■n'^ 

+ ± 2 ^ ^ + 3^r-2n2 + 


(4.4) 


The values of the sums over frequencies is given for various grids in the 
first Brillouin zone in Table 3. 

The harmonic approximation for the excess specific heat is 


Harm 


Nk 


= 1 


(4.5) 


B 


We can also calculate in the harmonic approximation root-mean-square 

H 27 

deviation of particles from their lattice sites, Y • It is given by 


qj “j 


j ^j(q) 

coth ) 

m 2N UJ,(q) 2 k^T 


(4.6) 


where d is the near neighbor distance. In the classical limit this becomes 


(y“)" = <//d^ ^1/31 


(jUq 2 

T 2rr N 5 ^uTTq) ^ 
qj j 


(4.7) 


This result is tabulated for varying numbers of points in the Brillouin zone 
in Table 5. 


An important point to note is that the small q behavior of cu.(q) is 

w 
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ujj_(q) “ /q 




(4.8) 


Thus, the harmonic approximation predicts a logarithmic divergence^ as the 
thermodynamic limit is approached, of due to the small q behavior of 
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(JUoCq.) • This: divergence has been seen in molecular dynamics studies 
d 

of crystalline systems of hard disks ranging in size from 100 to approximately 
7,000. We have no reason to believe that it will not occur in our system. 
This means that Lindemann' s ratio is not independent of the size of the system 
in two dimensions. We give a more detailed discussion of these matters in 
Section VI . 

Our solid phase MC results are tabulated in Table 4, together with some 
of the lattice dynamics predictions. Our MC results are well parameterized 
by adding a small ant^armonic correction to the harmonic excess internal 
energy (4.1) 


Harm 


ar”^ + br”^ 


(4.9) 


The anharmonic contribution is the first part of an expansion in powers of 
T or inverse powers of P. The values of a and b which we find are 5.0 and 
560. The internal energy is' displayed in Figure 1. 

Given our simple result for the excess internal energy the other thermo- 
dynamic functions follow directly. We assume that in the infinite P limit the 
harmonic free energy becomes the exact free energy. Hence, we may obtain the 
free energy for finite F by integrating the anliarmonic correction to obtain 


^tot harm 
F F 

NiTT Nk^t 
B B 


+ J (aP'"^+ bP'”^) 


03 



_harm 

£ 

NbT 



(4.10) 


The free energy is plotted in flgvire 2. By differentiating the excess internal 
energy according to Eq. (3,31) we obtain the excess specific heat of the solid 
phase, 


ex 

^ = 1 + 2aP ^ + 3bt . (4.11) 

* B 

Figure 3 shows this function and also our result for from direct Monte Carlo 


computation. 
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V. THE PHASE TRANSITION 

Using the free energies calculated in the preceeding sectLons we searched 
for a phase transition between the liquid and solid phases. Looking 
at table 2 we see that for the 100 particle system the liquid free energy is 
lower than the solid free energy for F below about 130. However, this crossover 
point is extremely sensitive to the free energies. An error of only 0.04% 
in the total free energy or 0.7% in the excess thermal free energy would shift 
the melting point F by 15. Doing a double tangent construction to determine 
the width of the two phase region shows that the melting and freezing points 
are only separated by about 0,1 in F. 

Since the free energies lie so close it is worthwhile to seek confirmation 
of the location of the phase transition by looking at the behavior of systems 
which were started from a perfect crystal and allowed to melt. The 100 parti- 
cle system started at F = 130 achieved an equilibrium value of the root mean 
square deviation of .138. This system was allowed to age for 13,000 moves per 
particle. The F = 120 particle system melted after about 2000 moves per 
particle. In our N = 64 particle runs the same behavior was observed with the 
r ss 130 system attaining an equilibrium root mean square deviation of 0.153. 
Hence, the phase transition probably lies below a value of F of 130 and above 
a F of 120. Since metastability appears to be much less of a problem with 
softer potentials in two and three dimensions and the free energies are 
difficult and expensive to compute for such potentials , monitoring the stability 
of the crystal lattice is a sensible alternative for soft-potentials. The 
quantity r is essentially a measure of the ratio of the potential and kinetic 
energies of the plasma. We therefore see that the system crystallizes when 
the potential energy is approximately one hundred times larger than the kinetic 
energy. This, to us, somewhat surprising result is very similar to that found 


in three dimensions. 
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Hockney and Brown^^ found a second order phase transition at F = 95 + 2 
in a molecular dynamics simulation involving 10,000 particles. We have made 
several careful runs in the neighborhood of F = 95 but have found no signs of 
any anomaly in either the free energy or the internal energy. In figure 3 
and table 2 we have presented the results of a direct calculation of the speci- 
fic heat at constant area based on the fluctuations in internal energy: 


<( 2 V(? ) ) > - < 2 V(r. ,))^1 (5.1) 

i<3 ^ i<J 

Calculation of the specific heat in this fashion is inevitably noisy but our 

11 

results are clearly incompatible with the results of Hockney and Brown . They 
performed their calculation by starting at a low temperature, .or high value of 
F, with a crystal with several grain boundaries and then increasing the tempera- 
ture in steps. We interpret the discrepancy with our results as showing that they 
did not give their system time to achieve equilibrium at the various tempera- 
tures. It is hard to blame the discrepancy in the difference in numbers of 
particles; we do not see how a system which melts at a value of F of 120 for 
100 particles could remain stable at a F of 100 for 10,000 particles. As the 
next section will show, traditional indicators of crystalline order, such 

as the mean square deviation of particles from their lattice sites, increase with 
the number of particles for fixed F. We know of no quantity which indicates 
increased order as the number of particles increases, in two-dimensional systems. 

The order of the phase transition is an important question. A recent 
discussion by Kosterlitz and Thouless argues that melting in two dimensions 
for short range forces is caused by the appearance of free dislocations as a 
result of the breakup of pairs (or higher combinations) of dislocations with opposing 
Burger's vector. Their calculation results in an analytic specific heat at the transi 
tion but the approximations they make are of the type which could easily mask 




a weak singularity. A recent calculation by Nelson , yields an essential 
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singularity in the specific heat. Another calculation by Hole and Medeires 
argues for a first order phase transition with short range forces and gives a 
rationalisation for the second order phase transition observed by Hockney 
and Brown^^ . 

Young “ has performed a calculation, similar to that of Nelson , paying 

particular attention to the angular forces between dislocation pairs, and finds 

33 

qualitatively similar results. Most recently Halperin and Nelson have argued 
that two-dimensional melting occurs in two steps. They propose that at a low 
temperature the breakup of dislocation pairs leads to a transition to a "liquid 
cx^ystal" phase and at higher temperature the dissociation of disclination pairs 
yields an isotropic fluid. 

We thus find the theoi*etical situation to be less than completely clear, 
especially for long range forces. Our results are compatible with a first 
order phase transition but our total free energy curves cross with a difference 
in slope of 0.03%* It may be argued that fitting the equation of state data 
biases one toward a first order phase transition and we are unwilling to state 
an order for our phase transition. All we can say is that in the thermodynamic 
quantities we have calculated we see no indication of any divergences . 

It is our conviction that the way to proceed at this point is to attempt 
to use molecular dynamics to investigate the mechanism of melting and we are 
starting further work along these lines. 

VI. LINDBliVNN'S R.ATI0 ..VND 3ISS DEPENDENCE 

In table 4 we show values of the root mean square displacement of electrons 
from their original lattice sites. All of these values are for S4 or 100 particle 
systems with periodic boundary conditions. Examination of these quantities 
shows that melting occurs for a root mean squai'e displacement of .16 + .01 in 
terms of the near neighbor distance. This quantity is known as Lindemann's 
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ratio. However, as we will see, this statement is an over slaplilicatioii of the 

situation in two dimensional crystals . 

"34 35 

Over 40 years ago Peierls and Landau argued that there would he no true 
long range two-dimensional crystalline order. Peierls produced an argument 
based on the harmonic approximation and Landau used his theory of second order 
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phase transitions. Ten years ago these arguments were made rigorous by Mermin 
who proved that, for every k 0, the Fourier component of the density must vanish 
the thermodynamic limit. More precisely he showed that 


p- (0n(N))”- 

where ^ ^ . 

K ^ i=l ’ 


( 6 . 1 ) 

( 6 . 2 ) 


He also showed that 


:|u(r)|M .. 


const N) 


(6.3) 


Where u(R) is the deviation of the particle from the lattice point at R. 
This proof is valid for potentials i2J<r) for which 


0(r) - Xr“|7“0(r) 


(6.4) 


is integrable at r = * and positive and nonintegrable at r - 0, both for 

X =t 0 and some positive finite value of X. The 1/r potential does not meet 

the first criterion. Hence, the question as to whether the two-dimensional 

one component plasma can display long range crystalline order has not been 

34 35 

rigorously answered at present. We do, however, find the Landau- Peierls ’ 
argument very convincing. 

To investigate this question we have performed a series of Monte Carlo 
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calculations at a value of F equal to 200 and with the number of particles varying 
from 16 to 1024. The root mean square deviation is plotted against (^(N))^ 
tn figure 4. For comparison the result of a lattice dynamics calculation is 
shown on the same graph. Due to the extremely lengthy calculations required 
to achieve convergence of this quantity we were not able to obtain more than 

a lower bound for the cases of N = 512 and N = 1024. It is seen that the 

points from N = 16 to N = 256 are compatible with the 0 ti(N) behavior but 
do not definitely rule out the approach of the root mean square deviation 

a constant value. Ve plan more work, on this question in the near future. 

In figures 5 and 6 we have displayed plots of the distribution of particles 
about their lattice sites for N values of 144 and 1024. The 
increase in the mean square deviation can be clearly seen. The one particle 

distribution functions are displayed in figure 8 . 
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In his paper Mermin pointed out that two-dimensional crystals, while 
not possessing true long range translational order, may have long range orienta- 
tional order. If R a ■** * R + n(R) then in the harmonic 

approximation 


= <[r(Rt-aj_) - r(R)]’ Cr(ISl-a^) - r(R)]> - | a^_ 1 as 1 r-R' | - co (6.3) 


In figure 7 we show the results for A"' at F = 200 for various numbers of 
particles . It is seen that this quantity rapidly approaches a constant 
value independent of N. 

We believe that we have established that the lattice displacements in 

the one component plasma behave in a very similar way to those of the hard 
28 

disk system . There is a loss of translational long range order as the size 
of the system increases. On the other hand there appears to be long range orienta 
tional order. The full implications of this observation must await a more 


detailed investigation of larger systems. 
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VII. ONE AND TWO PARTICLE DISTRIBUTION FUNCTIONS 

In. figure 3 we have plotted the distribution function, n(r) , for particles 
about their lattice sites. The single particle distribution function is defined 
30 that n(r)dr is the probability of finding a particular particle within a 
volume element dr at a point removed from the lattice site by a displacement 
r. The logarithm of n(r) has been plotted as a function pi the distance from 
the lattice site squared. Hence, if the distribution of particles about their 
lattice sites were Gaussian, the points would fall on a straight line. The 
normalization has been arbitrarily chosen so that n(r) = 1 for the point with 
the smallest value of r. 

Two effects are illastrated in this plot. For the two plots with N = 64 

wo can see the effect of lowering from 300 to 200. In addition, the effect 

of changing N from 64 to 1024 for T = 200 is evident. In none of the cases is 

n(r) truly Gaussian. The result for F = 300 is approximately Gaussian but the 

results for F = 200, although more spread out than n(r) for F = 300, appear to 

be cutoff more steeply than a Gaussian distribution. Our results contrast with 
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the results of Young and Alder* who find that densely packed hard disks form 
a Gaussian distribution but that at lower density the distribution decays more 
slowly than a Gaussian distribution. Thus, the considerable differences between 
hard disks and 1/r particles seem to produce opposite deviations from Gaussian 
behavior as the melting transition is approached. 

The radial distribution function g(r) is defined by the equation, 

j(r, „(r )]d? d? (7.1) 

p ^N i<j 

We have calculated values for g(r) function for F = 36, F — 90, and F = 120, 
These are tabulated in table 6 and plotted in figures 9, 10, and 11. For F = 90 
the HNC result for g(r) is plotted for comparison with the Monte Carlo result. 
The perhaps surprising feature is that g(r) shows considerably more structure 
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than the corresponding values of g(r) given by Hansen “ for the three-dimen- 
sional case. This, however, may be a general feature of two-dimensional 

36 37 

simple fluids as both Fehder and Tsien and Valleau found the height of 
th« first peak in g(r) to fall between 3 and *i for two-dimensional Lennard- 
Jones fluids. The HNC result for g(r) shown for T ” 90 also shows correspondingly 
more structure than three-dimensional HNC results. 

Finally, we have also calculated the structure factor S(k) which is 
defined via 

(7.2) 

(7.3) 


SW =-< pjp_-> 


Where 


^ « ikT, 

p. = 1 e i 

^ i=l 


The k vectors used are those corresponding with the reciprocal lattice generated 

via the periodic boundary conditions associated with the N-particle basic 

^tQnte Carlo rectangle. To determine S(k) we directly used the definition 

(7.2) with (7.3). In Table 7 and figures 12 and 13 our results for S(k) are 

illustrated for =: 36 and T = 90. For = 90 we also present the result of 

6ur HNC calculation. Just as in the case of g(r) these show more structure 

12 

than the three-dimensional results of Hansen . ^Ve found it difficult to get 
a good estimate of the height of S(k) at the first peak but it is clear that 
it will be larger than the value of 2.86 found at freezing in many three- 
dimensional systems*^' . 

The Debye-Hlickel result for S(h)^ 




(7.4) 


provides the correct low k behavior for S(k). However, for our smallest values 
of k the Monte Carlo values of 3(k) have already risen abwe the Debye-Huckel 




27 


VIII. CONCLUSIONS 

In this paper we have presented the results of a Monte Carlo calculation 
of th;‘ properties of particles interacting via the 1/r potential in two 
dimensions. In particular, we have emphasized the nature of the ordered 
phase and attempted to show how the dimensionality has influenced the nature 
of the order. The second main point is the phase transition itself. Much 
work is currently underway in two-dimensional melting and much remains to be 
done. We hope that experimentalists will soon observe an ordered phase of two- 
dimensional electrons. This should be possible with lower temperature experi- 
ments as higher values of can be achieved at lower densities. In addition, 

as figure 14 shows, quantum effects, as measured by the ratio of r and the 

s 

de Broglie thermal wavelength, will become less important in the neighborhood 
of the phase transition if it is observed at lower temperatures. 
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Figure Captions 
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Fig. 2 


Fig . 3 


Fig. 4 
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The excess thermal internal energy as a function of T. The 
points are the results of Monte Carlo calculations with the 
triangles representing a 16 particle system, the circles a 36 
particle system, the squares a €4 particle system and the 
diamonds a 100 particle system. 

The excess thermal free energy as a function of F- The circles 
give the results of a calculation based on the hypemetted 
chain integral equation, the triangles show the results of the 
Monte Carlo liquid state calculation, the diamonds show the 
predictions of lattice dynamics in the harmonic approximation, 
and the squares show the results of our Monte Carlo solid state 
Calculations. 

The excess specific heat as a function of F. The points give 
the result of direct Monte Carlo calculation based on equation 
(5.1). The solid line shows the result of calculating the 
specific heat via equations (3.33) and (4.11) which were obtained 
from fitting the Monte Carlo results for the internal energy. 

The root mean square deviations of particles from their lattice 
sites as a function of the square root of the logarithm of the 
number of particles. The root mean square deviation is measured 
in terms of the near neighbor distance, d. The solid circles 

I 

are the results of Monte Carlo calculations for 16, 36, 64, 100, 
144, 256, 576, and 1024 particles. The point for 1024 particles 
represents only a lower bound. The circles represent the predic- 
tion of lattice dynamics in the harmonic approximation. The 

squares represent the results of a molecular dynamics calculation 
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for hard disks by Young and Alder . The hard disk calculations 



Fig. 5 


Fig. 6 


Fig. 7 


Fig. 3 


have V/V = 1.20 where V is the close packed volume, 
o o 

This figure shows the superimposed positions of the particles 
for 100 different Monte Carlo configurations, each separated 
by 4 passes. A pass Is defined as one attempted move per 
particle. In this simulation F =» 200 and N = 144. Roughly 
speaking, this represents a short time picture of the crystal. 

This figure shows the superimposed positions of the particles 
for 100 different Monte Carlo configurations, each separated 
by 4 passes. A pass is defined as one attempted move per 
particle. In this case F » 300 and N » 1024. If this figure is 
compared with figure 5, the additional effect of the long wave- 
length phonons for 1024 particles may be seen. To see the long 
wavelength osclllatldns in particle positions look down the 
rows from a vantage point almost in the plane of the paper. 

The angular correlation tT t defined in equation (S.3), as a 

function of the distance, R, between pairs of particles for 
2 

F 3 200. Both ^ and R are measured in terms of the near neigh- 
bor distance. The results are for Monte Carlo calculations with 

2 

144, 256, 576, and 1024 particles. In each case ^ differs only 
very slightly from 1, thus demonstrating the angular order 
observed in the crystal. 

The single particle distribution function as a function of r where 
r is measured in units of the near neighbor distance. The triangles 
are for 1024 particles with F = 200, the squares are for 64 
particles with F =200, and the circles are for 64 particles 
F = 300. The normalization has been chosen so that n(r) = 1 


for the first point. 


Plff. 9 


n«. 10 


Fig, 11 


Fig* 12 


Fig. 13 


Fig. 14 


The radial distribution tunction, g(r) , as function of intBr~ 

particla distance. The points were produced via Monte Carlo 

calculation with 100 particles and are for T =* 36. The distance 

1 

is measured in terms of the ion sphere radius, a ^ l/(Tm) , 
where n is the density. 

The radial distribution function, g(r) , for T «=» 90, The triangles 
are the results of Monte Carlo simulation with 100 particles and 
the circles were calculated by solving the hypernettad chain 
equation. Distances are measured in terms of the ion sphere 
radius. 

The radial distribution function, g(r) , for T =* 120. The points 
are the result of Monte Carlo simulation with 100 particles. 
Distance is measured in terms of the ion sphere radius. 

The structure factor, S(k), for F » 36. The points resulted 
from Monte Carlo calculation with 100 particles. The wavenumber, 
k, is measured in terms of the inverse ion sphere radius. 

The structure factor', S(k) , for F = 90. The triangles are the 
result of Monte Carlo calculation while the circles were calcu- 
lated via solution of the hypernetted chain integral equation. 

The Wavenumber, k, is measured in terms of the inverse ion sphere 
radius. 

.4 density versus temperature plot which shows the location of the 
predicted phase transition. Only the region of temperature greater 
than about 1®K and density less than approximately 2 x lO^cra""" 
have been e^cplored experimentally. The line where the de Broglie 
thermal wavelength is one quarter of the near neighbor distance has 
been included to give an indication of the region where quantum 
mechanical effects become important. 


TABLE I 


Flttlog parameters for the liquid phase results 

c 3 >9.290414 



MC 

HNC 

®1 

-1.106103 

-1.102071 


.765873 

.799066 

^3 

.775443 

.951230 


.261904 

.201743 

^5 

-1.202048 

-1.593872 

*0 : 

.957986 

.131187 

V 

- 

- .232854 

*8 

- 

.536553 


Solid phase parameters: a s 4.9S6, b s 561.1 



TABLE II 


N 

r 

th , 

U /NkgT 

F^®^/NkgT 

_ex 

64 

1 

.32 

.61 

.16 

64 

2 

.43 

- 1.46 

.26 

64 

5 

.61 

- 4.30 

.40 

64 

10 

.75 

- 9.36 

.52 

64 

20 

.91 

- 19.84 

.66 

64 

30 

1.01 

- 30.52 

.76 

64 

40 

1.08 

- 41.28 

.82 

64 

50 

1.14 

- 52.09 

.39 

100 

60 

1.19 

- 62.93 

.94 

64 

70 

1.23 

- 73.81 

1.00 

100 

80 

1.25 

- 84.70 

.37 

64 

90 

1.29 

- 95.61 

1.16 

100 

100 

1.33 

-106.54 

1.06 

100 

110 

1.34 

-117.47 

1.36 

100 

120 

1.32 

-128.42 

1.37 

The 

thermodynamic functions 

of the liquid 

phase. The internal energy 

: had 

the ideal gas and static 

lattice contributions subtracted 

from it and 


the specific heat has had the ideal gas contribution subtracted. The number 
of particles, N, is the number of particles in the Monte Carlo run. 



TABLE III 


N 

2 *> 
o 



<^Oo/uj )> 
o 

4 

6.63 

3.13 

1.67 

-.599 

16 

3.S8 

3.39 

1.63 

-.677 

36 

10.38 

3.49 

1.63 

-.693 

64 

11.43 

3.53 

1.63 

-.698 ' 

100 

12.25 

3.56 

1.63 

-.700 

144 

12.91 

3.38 

1.63 

-.702 

196 

13.47 

3.59 

1.63 

-.702 

The results 

of our lattice 

dynamics ( 

calculation are presented 

in the 

following table. 

N represents 

the number 

of points summed over in 

the 


Brillouin zone. The an^lar prackata deilote the average of the various 

2 2 3 

functions of frequency. The characteristic frequency tu“ =* a“/ma . 
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TABLE IV 


N 

r 

th , 

U /NkgT 

F*°VNkgT 

ex 

C^/NkgT 

2 2 
<r /d > 

100 

130 

1.06 

-139.37 

1.25 

.16 

100 

140 

1.06 

-150.35 

1.15 

.15 

100 

150 

1.03 

-161.34 

.98 

.14 

100 

170 

1.04 

-183.33 

1.13 

.13 

100 

190 

1.04 

-205.34 

1.09 

.12 

64 

220 

1.03 

-238.37 

1.10 

.11 

64 

250 

1.03 

-271.42 

1.03 

.099 

64 

300 

1.03 

-326.54 

1.08 

.092 

The 

thermodynamic 

functions of the 

solid phase. 

The internal energy- 


had 

the ideal gas 

and static lattice 

contributions 

subtracted from it 

and 


the specific heat has had the ideal gas contribution subtracted. The number 

of particles, N, is the number of particles in the Monte Carlo run. The 
2 2 

column headed <r /d > gives the root mean square deviation of particles from 


their lattice sites. 



TABLE V 




16 

.09 

.11 

36 

.10 

.12 

64 

.11 

.12 

100 

.12 

.13 

144 

.12 

.13 

256 

.13 

.14 

512 

.13 

.14 

1024 

>.13 

.15 


This table gives the root meah square deviation i£ particles from 
their lattice sites for r = 200 and various numbers of particles in the 
basic MC rectangle* The predictions of lattice dynamics in the harmonic 
approximation are also presented. 



TABLE VI 


Monte Carlo 


r 

8(r; r =36) 

g(r; r =90) 

r 

g(ri r =36) 

sir; T>90) 

1.06 

.01 

.00 

4.31 

.89 

.71 

1.11 

.02 

.00 

4.37 

.89 

.71 

1.17 

.06 

.00 

4.43 

.39 

.71 

1.23 

.13 

.00 

4.48 

.39 

.72 

1.28 

.25 

.01 

4.54 

.91 

.76 

1.34 

.44 

.05 

4.60 

.93 

.79 

1.40 

.66 

.16 

4.66 

.95 

.84 

1.46 

.94 

.37 

4.71 

.95 

.89 

1.31 

1.22 

.73 

4.77 

.98 

.94 

1.57 

1.46 

1.21 

4.83 

1.00 

1.00 

1.63 

1.68 

1.72 

4.39 

1.02 

1.06 

1.69 

« 1.81 

2.13 

4.94 

1.03 

1.11 

1.74 

1.37 

2.49 

5.00 

1.05 

1.15 

1.80 

1.86 

2.63 

5.06 

1.05 

1.18 

1.86 

1.81 

2.55 

5.11 

1.06 

1.21 

1.91 

1.71 

2.33 

5.17 

1.07 

1.22 

1.97 

i -57 

2.06 

5.23 

1.06 

1.22 

2.03 

1.44 

1.75 

5.28 

1.06 

1.22 

2.09 

1.30 

1.45 

5.34 

1.06 

1.19 

2. l 4 

1.17 

1.20 

5.40 

1,05 

1.17 

2.20 

1.06 

.99 

5.46 

1.04 

1.14 

2.26 

.95 

.81 

5.31 

1.03 

1.10 

2.31 

.87 

.70 

5.57 

1.02 

1.06 

2.37 

.79 

.60 

5.63 

1.01 

1.03 

2.43 

.75 

.52 

5.69 

1.00 

.99 

2.48 

.70 

.48 

5.74 

.99 

.95 

2.54 

.67 

.45 

5.80 

.98 

.92 

2.60 

.68 

.44 

5.36 

.98 

.39 

2.66 

.66 

.44 

5.91 

.97 

.87 

2.71 

. 66 

.45 

5.97 

.97 

.35 

2.77 

.69 

.47 

6.03 

.96 

.35 

2,83 

.72 

.52 

6.09 

.96 

.35 

2.39 

.76 

.56 

6.14 

.96 

.85 

2.94 

.30 

.64 

6.20 

.97 

.87 

3.00 

.87 

.71 

6.26 

.97 

.89 

3.06 

.92 

.30 

6.31 

.97 

.91 

3.11 

.99 

.91 

6.37 

.98 

.93 

3.17 

1.04 

1.03 

6.43 

.98 

.95 

3.23 

1.09 

1.14 

6.48 

1.00 

.98 

3.23 

1.14 

1.25 

6.54 

1.00 

1.01 

3.34 

1.17 

1.35 

6,60 

1.01 

1.03 

3.40 

1.19 

1.44 

6.66 

1.01 

1.06 

3.46 

1.21 

1.48 

6.71 

1.02 

1,07 

3.51 

1.21 

1.49 

; 6.77 

1.03 

1,10 

3.57 

1.19 

1.48 

; 6.83 

1 . 02 

1,10 

3,63 

1.19 

1.45 

6.89 

1.02 

1.11 

3.69 

1.17 

1.39 

: 6,94 

1.02 

1,10 

3.74 

1.13 

1.31 

7.00 

1.02 

1.10 

3.80 

1.10 

1.22 

7.06 

1.02 

1.09 

3.36 

1.06 

1,14 

7.11 

1,01 

1.08 

3.91 

1.03 

1.06 

7.17 

1.01 

1.06 




i 1.1 

1 } T I r 

1 T T T ' 



TABLE VI 
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r 

g(r; r*36 

g(r; r-90) 

r 

g(r; r=36) 

g(r; r=90) 

3.97 

1.00 

.97 

7.23 

1.00 

1.05 

4.03 

.96 

.90 

7.28 

1.00 

1.03 

4.09 

.94 

.34 

7.34 

1.00 

1.01 

4.14 

.92 

.79 

7.40 

1.00 

.99 

4.20 

.91 

.75 

7.46 

.99 

.98 

4.26 

.90 

.72 

7.51 

1.00 

.97 



Monte Carlo 




Hypemetted Chain 

r 

g(r; ra 36) g(r; 

r*90) 

r 

g(r; r *120) 

r 

g(r; r =*90) 

7.57 

.99 

.94 

1.30 

.01 

1.17 

.01 

7.63 

.99 

.94 

1-45 

.20 

1.30 

.12 

7.69 

.99 

.93 

1.60 

1.31 

1.44 

.67 

7.74 

.99 

.92 

1-75 

2.72 

1.57 

1.53. 

7.80 

.99 

.92 

1.90 

2.68 

1.70 

2.20 

7.86 

.99 

.93 

2.06 

1.70 

1.84 

1.96 

7.91 

.99 

.93 

2.21 

.90 

1,97 

1.61 

7.97 

.99 

.94 

2.36 

.51 

2.10 

1.25 

8.03 

1.00 

.95 

2.51 

.36 

2.24 

.98 

8.09 

.99 

.97 

2.67 

.35 

2.37 

.80 

8.14 

.99 

.98 

2.32 

.44 

2.50 

.70 

3.20 

1.00 

.99 

2.97 

.63 

2,64 

.65 




3.12 

.94 

2.77 

.66 




3.28 

1.30 

2.90 

.72 




3.43 

1.56 

3.04 

.83 




3.58 

1.59 

3.17 

.99 




3.73 

1.41 

3.31 

1.14 




3.89 

1.12 

3.44 

1.24 




4.04 

.35 

3.57 

l .;25 




4.19 

. 66 

3.71 

1.19 




4.34 

.59 

3.34 

1.10 




4.49 

,64 

3.97 

1,02 




4.65 

.78 

4. 11 

.95 




4.80 

,99 

4.24 

.90 




4.95 

1.20 

4.37 

.88 




5.10 

1.32 

4.51 

.88 




5.26 

1.31 

4.64 

.91 




5.41 

1,20 

4.77 

.95 




5 .56 

1-06 

4.91 

1.00 




5.71 

.94 

5.04 

1.05 




5.87 

.83 

5.18 

1.07 




6.02 

.77 

5.31 

1.08 




6,17 

.78 

i 5.44 

1.06 




6.32 

.86 

5.58 

1.04 




6.48 

.98 

5.71 

1.01 




6 .63 

1.12 

5.84 

.98 




6.78 

1.19 

5.98 

.97 




6.93 

1.19 

6.11 

.96 




7.09 

1.12 

6.24 

.96 


TABLE VI 


(continued) 


Hypernetted Chain 


£ 

g(r; r=l20) 

r 

g(r; F=901 

7.24 

1.03 

6.38 

.97 

7.39 

.95 

6.51 

.98 

7.54 

.89 

6 . 64 

1.00 

7,69 

,87 

6.78 

1.02 

7.85 

.88 

6.91 

1.02 

S.OO 

.92 

7.05 

1.03 

8,15 

.99 

7.18 

1.02 



7.31 

1.02 



7.45 

1.00 



7.58 

.99 



7,71 

.99 



7.85 

.98 



7.98 

.98 



3.11 

,99 



3.25 

.99 



3.33 

1.00 


Monte Carlo and hjrpometted chain results for the radial distribution 
function, g(r) , for F = 36, 90, and 120. Distance is given in units of the ion 
sphere radius. 



TABLE VII 


Monte 

Carlo 


H 3 rpernetted 

Chain Equation 

k 

S ( k ; r «36) 

S(lc; r =90) 

k 

S(k; r =90) 

1.40 

.04 

.02 

0.75 

.01 

1.69 

.06 

.02 

0.90 

.01 

1.96 

.09 

.04 

1.06 

.01 

2.16 

.12 

.05 

1-22 

.02 

2.36 

,16 

.08 

1.37 

. 0 ’\ 

2.55 

.23 

.10 

1.53 

.03 

2.72 

.32 

.16 

1,69 

.04 

2.66 

.41 

.21 

1.85 

.05 

3.02 

.63 

.35 

2.00 

.06 

3.14 

.78 

.52 

2.16 

.09 

3.29 

1.11 

.33 

2.32 

.12 

3.43 

1.43 

1.46 

2.47 

.17 

3.52 

1.30 

2.04 

2.63 

.24 

3.66 

1.94 

3.37 

2.79 

.34 

3.78 

2.08 

3.36 

2.95 

.51 

3.88 

2.05 

2.39 

3.10 

.78 

3.99 

1.98 

2.41 

3.26 

1.18 

4.08 

1.81 

1.90 

3.42 

1.69 

4.22 

1.49 

1.38 

3.57 

2.14 

4.29 

1.38 

1.24 

3.73 

2.26 

4.38 

1.27 

.95 

3.39 

2.07 

4.50 

1.10 

.87 

4.04 

1.75 

4.59 

1.06 

.75 

4.20 

1-45 

4.66 

.97 

.69 

4.36 

1.23 

4.76 

.93 

.67 

4.52 

1.06 

4.34 

.91 

.64 

4-67 

.94 

4.98 

.31 

.56 

4.83 

.85 

5.03 

.81 

.60 

4.99 

.79 

5.15 

.78 

.53 

5,14 

.75 

5.27 

.76 

.57 

5.30 

.73 

5.36 

.75 

.53 

5.46 

.71 

5.75 

.74 

.54 

3.62 

.71 

5.58 

.74 

.57 

5.77 

.72 

5.74 

.76 

.62 

5.93 

.73 

5.91 

.76 

,63 

6.09 

.76 

6.10 

.30 

.71 

6.24 

.80 

6.36 

.87 

.82 

6.40 

.84 

6.77 

.97 

1.06 

6.56 

.90 




6,72 

.96 




6.87 

1.03 




7.03 

1.10 




7.19 

1.15 




7.34 

1.19 




7.50 

1,20 

Mcnte 

Carlo and hypernetted chain results 

for the structure 

factor, S(k), 

r a 36, 90, 

and 120. 

Wavenumber is measured in 

terms of the inverse ion sphere 


rat^lua 
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0.16 


0.14 


0.12 


0 .10 
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"O 
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